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Abstract. We consider a non- interacting Fermi gas in a combined harmonic and 
periodic potential. We calculate the energy spectrum and simulate the motion of the 
gas after sudden replacement of the trap center. For different parameter regimes, the 
system presents dipole oscillations, damped oscillations around the replaced center, 
and localization. The behaviour is explained by the change of the energy spectrum 
from linear to quadratic. 

PACS numbers: 03.75.Ss, 03.75.Lm 
1. Introduction 

In solid state physics, electrons in an atomic lattice are modelled by fermions satisfying 
the Schrodinger equation in a periodic potential. Alkali atoms in a periodic potential 
created by light have recently provided a new physical realization of the same 
mathematical model [H 12 El IH El E]- Experimentally, optical lattices are highly 
controllable and open a wide range of parameters for study, also the dynamics. For 
bosonic atoms, quantum transport such as Bloch oscillations have been observed [2]. 
The first such experiments on fermions in optical lattices have been realized very recently 

13 il El El- 

Unlike electrons in bulk solids, the atomic gases feel in addition to the lattice an 
overall confining potential, usually of harmonic shape. This harmonic trapping can be 
chosen weak and the system considered nearly homogeneous. On the other hand, it 
can be made significant in order to induce interesting quantum transport experiments: 
shift of the trap center may cause oscillatory motion of the gas 0111 El- Single-particle 
states and spectrum in a combined harmonic and periodic potential have been studied 
in [71 Ej , and the effect of the potential on quantum transport has been analyzed in jJl Ej 
using semiclassical analysis. We present here exact numerical treatment which gives a 
complementary point of view for understanding the problem, and provides a method 
applicable for a wide range of parameter values. Since we consider non-interacting 
fermions, the only physical assumptions we need are that the individual atoms obey 
the Schrodinger equation and that the equilibrium of the many-particle system can be 
described by the grand canonical ensemble. The physical model is considered in more 
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detail in Section |21 Details of the numerical implementation are discussed in Section |21 
Finally, in Section |3] the results of the simulations are presented and discussed, and 
compared with the experiments. 



2. Physical model 

We consider a one-dimensional optical lattice combined with a three-dimensional 
quadratic potential. The quadratic potential is weak in the direction of the lattice (axial 
direction) and tight in the remaining two (radial) directions. The full three-dimensional 
Hamiltonian operator for a single particle is 

tt3D "' V72 i-"- / 22| 22, 2 2\|'-^ \ 
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where uja and ujr are the axial and radial frequencies, respectively, U is the lattice height, 
and A is the lattice wavelength. We shall consider the one-dimensional system in the 
axial direction first and deal with the radial degrees of freedom later. Choosing the 
units of energy and length to be hjJa and ^J h/muja, respectively, and using the notation 
A = U/2 and = 47r/A, the one-dimensional Hamiltonian can be written as a sum 

H = H^ + H^ = --4^7 + -x^ + Acos(kx) 
2 dx'' 2 

of the lattice potential = A cos{kx) and the oscillator Hamiltonian 
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Let {(pn) and (e„) be the eigenf unctions and eigenenergies of H, respectively. Using 
the annihilation and creation operators a„, a]^ associated to the eigenfunctions we may 
write the many-particle Hamiltonian in the fermionic Fock space as 



Initially the system is in the statistical equilibrium and can be described by the grand 
canonical ensemble 

$(0) = — ^ (1) 

Xre-^(^-''^) 

where n is the chemical potential, N = the number operator, and (3 = i-^. 

The time evolution of the system is given by 

$(t) = e^^''*/^$(0)e-^^'*/^ 

where H"^ is the Hamiltonian governing the dynamics of the system. Of course, H'^ 
must be different from the Hamiltonian H used to calculate the initial equilibrium state 
$(0), otherwise the system will be stable. In the experiment, the displacement d of the 
harmonic trap gives rise to the new Hamiltonian 

1 1 

H'^ = ---r^ + - + Acos{kx), 

2 dx^ 2 
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which starts the time evolution of the system. Displacement of the trap center in 
the absense of the lattice leads to the dipole oscillations characteristic for a harmonic 
potential. 

We want calculate the time evolution of the expectation value of the position 
operator. The task is greatly simplified by the fact that we can represent any many- 
particle state $ by a single-particle state p($) when calculating expectation values of 
single-particle observables. This reduction is discussed in more detail in the Appendix. 
The single-particle state p($(0)) corresponding to the grand canonical ensemble $(0) 
is a diagonal operator in the energy eigenbasis. The diagonal elements are given by the 
Fermi-Dirac distribution 

When the Hamiltonian operator is itself a single-particle operator, which is the case for 
non-interactig Fermi gas, the reduced state obeys the ordinary Schrodinger equation. 
Hence the final formula for the motion of the Fermi gas is simply 

{X){t) = Tr(e*^'*/V('^'(0))e-'^'*/^X), (2) 

where X is the position operator. 

As the final step we consider dimensional reduction from three dimensions to one. 
Even though we are interested in the motion of the Fermi gas only in one dimension, 
the physical system is naturally three-dimensional. However, we can simply trace out 
the radial degrees of freedom since they are independent of the observable which we are 
measuring. In practice, this is done by replacing the Fermi-Dirac distribution by the 
effective one-dimensional energy distribution 

oo 

feffie) = fi^+ik + l)hWr). 

k,l=0 

where the multiples of hur exhaust the energy spectra of the independent oscillators in 
the radial directions. 

3. Numerical implementation 

The time evolution of the center of mass of the non-interacting Fermi gas trapped in an 
optical lattice was calculated using the formula The first task was to diagonalize 
the Hamiltonian. Actually there were two Hamiltonians H and H'^ to be considered 
corresponding to the two positions of the trap. However, we made the simplifying 
assumption that the displacement of the trap was a multiple of the lattice wavelength 
(this is justified when the latter is much smaller than the former). Then by a symmetry 
argument the two sets of eigenfunctions were obtained from each other by a simple 
translation. According to the formula (PJ) the initial equilibrium state $(0) is diagonal 
in the eigenstate basis of the initial Hamiltonian H, and the diagonal elements are 
determined by the corresponding eigenvalues. The chemical potential fi was fitted to 
match the particle number. By a change of basis, the state matrix was then transformed 
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to the eigenstate basis of the displaced Hamiltonian H'^. The position operator X was 
also expressed in this basis. By definition, the Hamiltonian itself H'^ is diagonal in this 
basis, so that the time evolution was easy to calculate using formula (j2I). 

The main difficulty was to find a basis suitable for diagonalizing the Hamiltonian 
numerically. We used the Hermite functions 

where Hn is the nth Hermite polynomial. With this choice the harmonic oscillator 
Hamiltonian becomes diagonal with matrix elements 

i/l = (n + l/2)5„^ 

and even the matrix elements of the lattice part can be calculated analytically. To 
be exact, the matrix elements are recovered as = AIie{Inm) where 

ikx 



Inm = Inmik) := / (x)e* dx. 

J —oo 

We were able to derive explicit closed expressions for Inm{k) as finite polynomials of 
k. However, they are rather cumbersome and not very useful in practice. In numerical 
calculations the integrals are most easily obtained from the formula 

Inm = e-'^'\-^- (3) 
where hnm satisfies the recursion equation 

)/n (4) 

with the initial conditions /iqo = 1 and hnm = if either of the indices is negative. We 
proved formulas (jH)) and Q using generating function methods. The recursion formula 
is efficient but numerically unstable. Therefore high precision arithmetics had to be 
used. 

Since Fourier transform JF is a unitary transform on the space of square integrable 
functions and since Hermite functions are eigenf unctions of JF, J^ipn = i'^ipn, we can also 
use the above formulas to calculate the convolutions 

i^n{x)iJrn{x - d) dx = i'^^'^ I nm{d) . 

Apart from being interesting in their own right, these convolutions give the matrix 
elements of the change of basis from (ipnix)) to the translated basis {ipn{x — d)), which 
is exactly what was needed to transform the state matrix from the eigenstate basis of 
H to that of H'^. Hence formulas Q and ^ were in double use in the implementation. 

For some values of the parameters it was beneficial to use the scaled Hermite 
functions (-\/A?/'„(Ax)) with A > 1 in the diagonalization to improve resolution. Even 
in this case the matrix elements needed are easily derived from the above formulas. 
In the simulations we used less than 3000 basis vectors and calculated typically a few 
hundred eigenstates. This was probably not enough for all the cases but at this point 
our computer resources became a limiting factor. 
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Figure 2. Energy spectrum of a single atom in the combined one-dimensional potential 
at the lattice height s = (grey line), s = 3 (black line), and s = 8 (dashed line). The 
energy unit is the energy quantum of the axial harmonic oscillator. The ground state 
energies are set to zero. A transition from a linear spectrum (harmonic oscillator) to a 
quadratic one occurs when the lattice is imposed. The quadratic part of the spectrum 
corresponds to localized states. Non-localized states appear again for higher energies, 
see the black line which becomes linear for n > 450. 



4. Results 



We used the above methods to simulate the motion of the free Fermi gas of potassium 
atoms in an optical lattice confined in a magnetic potential and compared them to the 
experimental results described in El El ■ The average temperature of the gas in the 
experiment described in [H| was T = 90 nK, the lattice wavelength A = 863 nm, the 
axial and radial frequencies Ua = 27t x 24 s~^ and u;^ = 27r x 275 s~^, respectively, and 
the average atom number 25000. Recall that k = Att/X. The depth U = 2A = sE^ 
of the lattice given in the units of recoil energy = h? /2m\^ was varied. Initially 
the trap was replaced by 15 /im to excite oscillations. The results of the simulation are 
shown in Figure ^ The energies of the lowest lying single-particle eigenstates are also 
plotted against the ordinal number in Figure El to give some qualitative insight to the 
system. 
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In the absence of the lattice potential (s = 0) the system is reduced to a harmonic 
oscillator and the average position oscillates freely showing the characteristic dipole 
oscillation in a harmonic potential. None of the eigenstates is localized and the 
corresponding energy spectrum is linear. The energy quantum or the difference between 
the consecutive eigenvalues gives the frequency of the oscillations. 

When the height of the lattice is increased to s = 3, both localized and non- 
localized eigenstates play a role. The oscillations of the system are damped, which can 
be understood as follows. According to the formula 0, the center of mass of the cloud is 
a sum periodic terms whose frequencies are given by the differences (e„ — Em) / ^ between 
the eigenvalues e„ of the dynamical Hamiltonian. When the eigenvalues are not evenly 
spaced, destructive interference causes dephasing. 

Increasing the lattice height also increases the offset of the dipolar oscillations since 
more and more particles become localized. When the height of the lattice is increased to 
s = 8, localized eigenstates dominate. The system is almost frozen to its initial position. 
The energy spectrum becomes almost quadratic. The spectrum corresponds to states 
that are localized in individual lattice potential sites. Neighbouring sites are shifted 
in energy according to the superimposed harmonic potential, therefore the quadratic 
form of the spectrum (difference in energy of two localized particles depends essentially 
only on their locations in the lattice). Actually each energy is doubly degenerate (the 
symmetric and the antisymmetric states formed of lattice sites at opposite sides of 
the center) which can be resolved by a closer look at the spectrum. Therefore, to be 
accurate, localized particles correspond to a superpositions of these doubly degenerate 
states. Representative examples of localized and non-localized eigenstates are shown in 
Figure 01 

It is noteworthy that in the intermediary region lowest-lying states are not localized. 
This is indicated by the linear parts in the spectra near origin. Hence the atoms do 
not tend to localize when the system is cooled down to ultra-low temperatures, on the 
contrary. Note that also higher energy states may be unlocalized, i.e. have a linear 
spectrum, as is shown for one set of parameter values in Figure El This is associated 
with the existence of higher bands of the optical lattice potential. Such states may affect 
the localization and oscillation behaviour of the gas significantly in higher temperatures. 

The results are in very good qualitative and quantitative agreement both with the 
experiment and the semiclassical approximation presented in We compared our 
simulations also to two earlier experiments jBli^ by Modugno et al. In these cases the 
experimentalists observed oscillations with considerably larger amplitude and higher 
frequency than what was predicted by the simulations. We believe that this was simply 
because the physical parameters used in these experiments were much more demanding 
for our numerical method and would have required more eigenstates to be calculated 
than was currently possible. The larger the displacement and the particle number, and 
the higher the temperature, the more eigenstates are needed for reliable results. The 
computer resources at our disposal were rather modest so that we have hardly pushed 
the method to its limits. 
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Figure 3. The two cigcnstates (a) and (c) of the lowest energy associated to the 
potential (e) with parameters A: = 3 and A = 10 form a quasidegenerate pair. In 
a superposition of such states, the particle is localized in one of the wells of the 
potential (e). Relating to Figure 2, such states correspond to a doubly degenerate 
pair of eigenvalues in a quadratic spectrum. The eigenstates (b) and (d) associated to 
the potential (f) with parameters k = 12 and A = 30 are not localized. They resemble 
the eigenstates of the harmonic oscillator, modulated by the lattice potential. They 
correspond to a linear spectrum and their energies are not degenerate. 



5. Conclusions 

Wc have analyzed the energy spectrum and motion of a Fermi gas in a combined 
potential. The harmonic potential alone has a linear spectrum and shows dipole 
oscillations of the gas when the trap center is displaced. Introducing the periodic 
potential transforms the spectrum gradually into quadratic one and the gas becomes 
localized. Intermediate cases show damped oscillations around the displaced center 
position. Our approach is complementary to the semiclassical analysis presented in 
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5 S'lid is in good qualitative agreement with the experiments. For the parameter 
values which are not too challenging from the numerical point of view the results agree 
also quantitatively with both the experiments and the semiclassical analysis. 

It is important to note that the quantum transport in the combined harmonic and 
periodic potential is qualitatively rather different from the Bloch oscillation or Wannier- 
Stark ladder behaviour in tilted periodic potentials. The two cases approach each other 
only when the gas moves in a region where the harmonic trap can be approximated by a 
linear potential. This requires a much shallower trap than that used in the experiments 
13 in E] • Understanding the behaviour of the non- interacting Fermi gas sets the basis 
for investigating effects of interactions IHUDI and eventually superfluidity [HI [121 • 
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Appendix 

In this appendix we formulate the reduction from many-particle to single-particle state 
in the language of symplectic geometry. This formulation may be helpful in using the 
reduction in a more systematic way and perhaps for more general purposes. 

Let us first recall some concepts of classical mechanics using a point-like particle 
moving in three-dimensional space as an example. The configuration space of the particle 
is just C = whereas the coordinates of phase space M ~ x include the 
momentum p as well as the position x. In mathematical terms, the phase space is the 
cotangent space T*C of the configuration space and it comes equipped with a canonical 
non-degenerate two-form 



which is known as the symplectic form while the pair (M, uj) is referred to as a symplectic 
manifold. The symmetries of the phase space are the diffeomorphisms of M preserving uj. 
Accordingly, infinitesimal symmetries are the vector fields A on M whose Lie derivative 
satisfy L^i^ = 0. Then each infinitesimal symmetry generates a one-parameter group of 
symmetries at least locally, where the correspondence is given by the flow of the vector 
field. 

Under some mild conditions we can associate a classical observable or a function 
Ha on the phase space M to an infinitesimal symmetry A. The observable is defined 
by the equation 



3 




k=l 



uj{A,B) = dHA{B) 
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which is required to be vahd for all vector fields B on M. In particular, since the time- 
evolution of a classical mechanical system preserves the symplectic form, it is generated 
by an observable, which is know as the Hamiltonian. In the example above it is clear 
that both translations and rotations arc symmetries of the phase space. An infinitesimal 
translation is just a three- vector t G M^, and the corresponding observable is 

Ht{^,p) = p-t. 

Similarly, an infinitesimal rotation can be given by a three-vector r as r ■ R where 
the components of R are the infinitesimal rotations about the coordinate axes. The 
associated classical observable is then 

i7r(x,p) = (x X p) • r. 

The general pattern emerging is that given a Lie algebra Q of infinitesimal symmetries, 
the collection of associated observables is most conveniently expressed as a map p from 
the phase space M to the dual Q*, that is, the space of linear maps g — > M. In this way 
we associated the momentum p to the translations and the impulse moment x x p to 
the rotations. Therefore the map p : M — > g* defined by the equation 

{p{x),A) = Ha{x) 

for all A G g and a; G M is known as the moment map. The Lie group G corresponding 
to the Lie algebra g acts naturally on both g and g*, and the moment map is equivariant 
in the sense that p{gx) = gp{x)foT all g in G. 

Returning to quantum mechanics, let 7i be a Hilbert space, on which a compact 
group G acts unitarily. This is the general framework for a quantum mechanical system 
with the symmetry group G. Now the projective space 

F{n) = {ven \ \\v\\ = i}/u{i) 

of pure states is canonically a symplectic manifold, on which G acts preserving the 
symplectic form. Hence from the point of view of symplectic geometry, all that we have 
said previously applies. 

Let g be the Lie algebra of G. The observable associated to an infinitesimal 
symmetry ^4 G g is calculated to be the expectation value Ha{v) — {v\A\v) where 
A stands for the representation of A in H. There is also the moment map 

p : Fin) ^ g* 

from the projective space to the dual of g defined by the equation 
Ha{v) = {p{v),A). 

Hence the expectation value of A in the vector state v is given by (A) = (p(f ), A). We 
can extend the moment map to any state $ by linearity. Then the moment map image 
p($) is a kind of reduced state which can be used to calculate expectation values for all 
observables in g. The point is that p($) may be a much simpler object than the state 
$ we started with. Furthermore, the equivariance of the moment map means that 

p{U^U-^) = Up{^)U-^ 
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for every element U of G. In particular, if the Hamiltonian H itself belongs to g, then 
the time evolution of the reduced state is given by the usual Schrodinger equation. 

In the application we have in mind G is the group of unitary transformations 
of the single-particle Hilbert space V and q may be identified with the single-particle 
observables or the Hermitian operators on V. Using the pairing [B, A) = Tt{BA) we 
may identify g with its dual. Let Ti. be the fermionic Fock space associated to V. Then 
G acts naturally on 7i and by the general theory 



for any many-particle state $ and single-particle observable A. It is easy to check that 
p($) is a positive operator with trace 



since the number operator is indeed the representation of / G g in the fermionic Fock 
space. 

In conclusion, we can represent any many-particle state $ by a single-particle state 
when calculating expectation values of single-particle observables. It is easy to 
calculate explicitly in a given basis (vn) of V using equation lA.ll In terms of the 
annihilation and creation operators associated to this basis, the matrix elements of 
are 



is the Fermi-Dirac distribution which is usually derived using the partition function 
approach. The preceeding discussion was not just for the sake of an appealing 
generalization but also to show that the reduction to the single-particle state formalism 
gives not only the energy distribution but a state matrix obeying the Schrodinger 
equation. 
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